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ABSTRACT 

We investigate the fragmentation criterion in massive self-gravitating discs. We present 
new analysis of the fragmentation conditions which we test by carrying out global 
three-dimensional numerical simulations. Whilst previous work has placed emphasis 
on the cooling timescale in units of the orbital timescale, /?, we find that at a given 
radius the surface mass density (i.e. disc mass and profile) and star mass also play a 
crucial role in determining whether a disc fragments or not as well as where in the 
disc fragments form. We find that for shallow surface mass density profiles (p < 2, 
where £ oc R~ p ), fragments form in the outer regions of the disc. However for steep 
surface mass density profiles (p > 2), fragments form in the inner regions of a disc. In 
addition, we also find that the critical value of the cooling timescale in units of the 
orbital timescale, /3 C rit) found in previous simulations is only applicable to certain disc 
surface mass density profiles and for particular disc radii and is not a general rule for 
all discs. We find an empirical fragmentation criteria between the cooling timescale 
in units of the orbital timescale, /?, the surface mass density, the star mass and the 
radius. 

Key words: accretion, accretion discs - protoplanetary discs - planets and satellites: 
formation - gravitation - instabilities - hydrodynamics 



1 INTRODUCTION 

The evolution of protoplanetary discs has been explored at 
great length in the past to understand the processes by 
which they form in the early Class phase, the accretion 
from the molecular cloud core onto the disc in the Class 
I phase, the mass and angular momentum transfer in discs 
once they have formed as well as in the early Class II stage of 
an isolated disc, through to the disc dispersal mechanisms. 
One such concept that has been considered is the impor- 
tance of the disc self-gravity, particularly in the earlier pe- 
riod of a disc's lifetime. It is during the early stages when 
it is massive enough to be self-gravitating that the concepts 
of angular momentum transport and fragmentation in these 
discs are important. This is a particularly significant aspect 
when considering whether gas giant planets could form via 
this method. Historically, planet formation by gravitational 
instability has not been thought likely since the planets that 
form in this way are predicted to do so far out in a disc 



(> O(100) AU e.g. 


Rafikov 2005 


Matzner & Levin 2005 


Rafikov 2009; Clarke 2009 Bolcy 2009, Stamatellos & Whit- 


worth|2008||Kratter, Murray-Clay & Youdin 2010 


1, whereas 



until recently, giant planets have been found only at small 
radii (< 10 AU). 

Recent advances in observations, such as the discoveries 
of planets at large radii (O(10 — 100) AU) from the parent 



star (Kalas et al. 2008 Marois et al. 2008) call for an al- 



ternative planet formation mechanism other than the stan- 
dard core accretion method to be considered. Furthermore, 
theoretical advances, such as the increased likelihood that 
planets can form by gravitational instability closer to the 



parent star in low metallicity environments (Meru & Bate 



2010) where core accretion finds it difficult (Kornet et al 



20051, call for planet formation by gravitational instability 
to be scrutinised in much more detail. 

There are two quantities that have historically been 
used to determine whether a disc is likely to fragment. The 
first is the stability parameter ( Toomre|1964 |, 



ttEG' 



(1) 



where c s is the sound speed in the disc, K ep is the epicyclic 
frequency, which for Keplerian discs is approximately equal 
to the angular frequency, £1, E is the surface mass density 
and G is the gravitational constant. Therefore, once the sur- 
face mass density and the rotation of the disc have been 
established, the stability is purely dependent on the disc 
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temperature. Toomre (19641 showed that for an infinitesi- 
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mally thin disc to fragment, the stability parameter must 
be less than a critical value, Q C rit ~ 1. 



Gammie (2001 1 showed that in addition to the stability 



criterion above, the disc must cool at a fast enough rate. Us- 
ing shearing sheet simulations, he showed that if the cooling 
timescale can be parametrized as 



P ^cool^, 



where 



^cool — V* 



/dUcoolV 

\ dt ) 



(2) 



(3) 



u is the specific internal energy and dw coo i/di is the to- 
tal cooling rate, then for fragmentation we require /3 < 
3, for a ratio of specific heats 7 = 2 (in two dimen- 
sions). Rice, Lodato & Armitage (20051 carried out three- 



dimensional simulations using a Smoothed Particle Hydro- 
dynamics (sph) code and showed that this cooling parame- 
ter is dependent on the equation of state. They showed that 
fragmentation can occur if /3 < /3 cr it where /3 cr it ~ 6 — 7 for 
discs with 7 = 5/3 and /3 cr it ~ 12—13 for discs with 7 = 7/5. 



Gammie (2001) and Rice et al. (20051 also showed that in 



a steady state disc where the dominant form of heating is 
that due to gravitational instabilities, since the gravitational 
stress in a disc can be linked to the cooling timescale in the 
disc using 



Qgi = 



(4) 



1 

9 7(7-1)0 

the rapid cooling required for fragmentation can also be 
interpreted as a maximum gravitational stress that a disc 
can support without fragmenting, which they show to be 
Q^GI.max « 0.06. 

The concept of a fast cooling needed for fragmentation 
is very clear from previous work. However, the value of the 
critical cooling timescale, /3 cr it (and therefore, by equation)!} 
QJGi.max), does not appear to be too clear cut: Rice et al. 



( |2003[ ) found that for a 0.1M Q disc with surface mass den- 
sity profile, E oc J? -7 / 4 , extending to a radius, R out = 25 AU 
around a IMq star, the disc fragments using /3 = 3 but not 
for ft = 5, whereas for a disc with mass Mdi sc = 0.25M Q , 
the disc fragments for /3 = 5. On the other hand, [Rice et al.] 
( 2005 1 suggest that the fragmentation boundary is indepen- 



dent of the disc to star mass ratio. Clarke, Ha rper-Clark| 
& Lodato (20071 showed that the critical value of /3 (be- 



low which fragmentation will occur if the stability criterion 
is met) may depend on the disc's thermal history: if the 
timescale on which the disc's cooling timescale is decreased 
is slower than the cooling timescale itself (i.e. a gradual 
decrease in 0) then the critical value may decrease by up 
to a factor of 2. More recently, |Cossins, Lodato fc C larke 
(20101 showed that the critical value varies with the tem- 



perature dependence of the cooling law. In addition, they 
carry out a simulation of a self-gravit ating disc with su rface 



Rice et al. 



(2005) who 
and 



mass density profile, E oc R~ 3 ^ 2 (c.f. 
used E oc i? _1 ), with ratio of specific heats, 7 = 5/3 
show that the critical value /3 cr it ~ 4. Using equation |4j this 
is equivalent to QfGi.max = 0.1 which brings the result of 
CKGi.max = 0.06 described above into question. Yet a number 
of papers have been produced which base their work on the 
concept of a single critical value of /3 (or equivalently, a max- 
imum gravitational stress value) (e.g. [C larke 2009, Rafikov 
2009[ |Cossins et al.||2010| |Kratter et al.|2010[ ). 



However, there is a more fundamental question that 
arises from previous numerical studies using a parametrised 



cooling method. In many simulations (e.g. Rice et aL]|2003 
|2005 Clarke et al.|2007 1 , the fragments all form in the outer 



parts of the discs. If the fragmentation criterion for a self- 
gravitating disc only depends on f3, then if fragments form, 
one would expect them to form at all radii since all radii have 
the same value of /3. This implies that the cooling timescale 
is not the only parameter on which fragmentation depends. 

In this paper, we investigate the criteria for fragmenta- 
tion in greater detail. In Section[2]we analytically investigate 
how fragmentation may be expected to depend on various 
disc parameters. We test this analytical theory by carry- 
ing out global three-dimensional simulations, the numerical 
setup of which we describe in Section [3] In Section [4] we 
make initial comparisons between our simulations and pre- 
vious simulations by Rice et al. ( 2005[ ) as well as discussing 
the implications of the disc setup. In Section [5] we present 
our simulations, the results of which we describe in Section 
We finally discuss and make conclusions in Sections [7] and 
respectively. 



2 ANALYTICAL VIEW 

As discussed in Section [TJ for a disc to fragment one criteria 
is that the Toomre stability parameter, Q < 1. Making the 
approximation that K cp ~ = \J GM+/R 3 and using H = 
Cs/Q, where H is the isothermal scale height of the disc, the 
Toomre stability criterion becomes a condition on the aspect 
ratio, H/R: 



(5) 



(6) 



H ttT,R 2 
~R ~ M* ' 

Approximating the surface m ass density a s E 
Mdisc /(t R 2 ), equation [H] becomes ( Gammie|2001 1 

H < Afdisc 
R ~ M* ' 

The surface mass density of a disc can also be written in 
the form of a power-law, E = E (i? / R) p , where E is the 
surface mass density at radius R , and p is a constant for 
any one disc. Substituting this into equation[5] the condition 
for fragmentation becomes: 

H_ < TVER 2 _ 1VE RP p (2-p) 

R~ M+ ~ GM+ ' 1 ' 

Equations [6] and [7| show the following: 

(i) Increasing the disc mass or decreasing the star mass 
is likely to promote fragmentation since a greater portion of 
the disc is likely to fulfil the above criteria. 

(ii) The surface mass density profile may play a part 
in the fragmentation of a disc. If p = 2, the disc is scale- 
free and each radius is equivalent to any other radius: the 
right hand side (RHS) of equation [7] is constant with ra- 
dius and the ratio of the cooling timescale to the orbital 
timescale, P, is also a constant. Therefore, if the disc settles 
into a quasi-steady state where the internal heating due to 
the gravitational instabilities balances the cooling, we ex- 
pect Q also to be constant with radius (i.e. the left hand 
side, LHS, of equation [7j H/R, is also a constant). Conse- 
quently, if p = 2, either the entire disc should settle into 
a quasi-steady state or the entire disc should fragment. We 



Fragmentation criteria in self- gravitating discs 3 



note that in a fragmenting disc, because the heating, cool- 
ing and fragmentation timescales are all proportional to the 
dynamical timescale in the disc, the fragmentation should 
occur first (in absolute terms) at small radii. 

For p < 2, the RHS of equation [7] increases with radius. 
Although, H/R is likely to increase with radius as well, since 
H/R will typically increase more slowly than the RHS of 
equation [7] this condition is more likely to be satisfied in 
the outer regions of a disc. Conversely, for p > 2, the RHS 
of equation [7] decreases with radius and hence the condition 
is more likely to be satisfied at small disc radii. 

(iii) For a p < 2 disc with a low enough ft such that 
it can fragment, the exact value of ft may determine how 
much of the disc satisfies condition [jj If the cooling in a 
disc is fast such that ft is small, the temperature and hence 
sound speed, c s will decrease more rapidly than in a disc 
where ft is higher. Consequently, since H oc c s , the aspect 
ratio will be lower at any particular radius and hence the 
disc is more likely to satisfy condition [7] for smaller ft. Since 
gravitational instability typically develops on a dynamical 
timescale, td yn oc I/O, oc R 3 ^ 2 , the instability will develop 
in the inner regions first and therefore fragmentation will 
first occur as close to the inner regions as possible where the 
fragmentation criteria are satisfied. Since more of the disc 
satisfies the above criteria for decreasing ft, the radius at 
which the first fragment forms will also decrease. We there- 
fore expect the fragmentation radius to move inwards with 
more efficient cooling. 

(iv) Crucially, equation [jj shows that the radius is im- 
portant and that for a shallow surface mass density profile 
(p < 2), there does not appear to be a limit for fragmen- 
tation (if the disc cools fast enough): provided the disc is 
large enough, condition [jj will be satisfied (since typically, 
an increase in H/R with radius will be much smaller than 
the increase in the RHS of equation [jj. 



3 NUMERICAL SETUP 

Our simulations are carried out using an SPH code originally 



and 0.2, respectively, which are fixed throughout the simu- 
lations. Furthermore, our work will begin with a comparison 



developed by Benz ( 1990 \ and further developed by Bate, 
|Bonnell fc Price| Hl995| T~and |Price fc Bate] |2007j. It is a 
Lagrangian hydrodynamics code, ideal for simulations that 
require a large range of densities to be followed, such as 
fragmentation scenarios. 

We include the heating effects in the disc due to work 
done on the gas and artificial viscosity to capture shocks. 
The cooling in the disc is taken into account using the cool- 
ing parameter, ft (equation [2JI, which cools the gas on a 
timescale given by equation (3] 

In order to model shocks, SPH requires artificial viscos- 
ity. We use a common form of artificial viscosity by |Mon-| 
aghan & Gingold (19831, which uses the parameters «sph 



and ftspH- A corollary of including artificial viscosity is that 
it adds shear viscosity and causes dissipation. If this vis- 
cosity is too large, the evolution of the disc may be driven 
artificially, while if it is too small, it will lead to inaccu- 
rate modelling of shocks (Bate 1995]) . As discussed in |Meru| 
& Bate (20101, various values of the SPH parameters have 
been tested and we find that a value of q:sph ~ 0.1 pro- 
vides a good compromise between these factors. Since typ- 
ically, ftspH ~ 2qsph, we choose osph and ftspn to be 0.1 



with Rice et al. ( 2005 1 and so we use the same values used 



by them. We use an adiabatic equation of state with ratio 
of specific heats, 7 = 5/3. 

3.1 Numerical effects on fragmentation results 



| Rice et aTj ( |2005| showed that for a disc with ratio of specific 
heats, 7 = 5/3, the critical value of the cooling timescale in 
units of the orbital timescale required for fragmentation is 
/Scrit ~ 6 — 7. The SPH code used for the simulations pre- 
sented in this paper differs in the way the smoothing length 



of the particles is set from that of Rice et al. (20051: whilst 



their code sets the smoothing length by approximately fixing 
the number of neighbours that each particle has to ~ 50, the 
current version uses a variable smoothing length which does 
not fix the number of neighbours but allows the smoothing 
length to be spatially adaptive whilst maintaining energy 
and entropy conservation ( Springel fc Hernquist|2002 Mon- 
aghan 20 02]), wit h our particular implementation described 
by iPrice & Batel (120071). 



All explicit hydrodynamical simulations must limit the 
timestep based on the Courant condition. Our SPH code also 
applies a force condition and a viscous timestep condition 



(see Monaghan 1992 for a review). In the simulations pre- 
sented here, since we apply an explicit cooling rate, it is 
important to ensure that the following timestep criteria is 
also met: 

,/3 



At s; C 



9: 



where C is a constant less than unity. We investigate the 
effects of varying the constant C on the critical cooling 
timescale /3 cr it by testing values of C = 0.3, 0.03 and 0.003. 
We find that this does not have a significant effect on the 
fragmentation results and so we use C — 0.3 for the simu- 
lations presented here. However, the timestep criterion may 
become more important for small ft or for particles at small 
radii. Therefore, for those simulations carried out with small 
values of ft (^ 3) or where fragmentation occurs at small 
radii (< 5 AU), the simulations have been repeated with 
C — 0.03 to confirm that this does not play a part in the 
results. 



4 BENCHMARKING SIMULATIONS 

Table [I] summarises the parameters and fragmentation re- 
sults of the simulations presented here. Each simulation was 
run either beyond the point at which the disc attained a 
steady state (for > 6 outer rotation periods, ORPs), or un- 
til it fragmented. 



The simulations presented by Rice et al. ( 2005 1 also 



used an SPH code. However, since the way the smoothing 
length is set in our code differs to the code used by |Rice| 
et al. (20051, and since it is uncertain as to whether their 



timestepping considered the cooling timescale, we simulate 
the same disc that |Rice et al.| |2005) simulated in order 
to initially find the critical cooling timescale in units of 
the orbital timescale, ftait- This is done by setting up a 
1 M Q star with a 0.1 M disc made of 250,000 SPH parti- 
cles, spanning 0.25 ^ R ^ 25AU. The initial surface mass 
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Simulation 
name 


P 


V 


1 


Qmin 


Initial Q profile 


Fragments? 


Bcnchmarkl 


6 


1 


0.5 


2 


Decreasing Q 


n 


Benchmark2 


5.5 


1 


0.5 


2 


Decreasing Q 


y 


Benchmark3 


5.6 


1 


0.5 


2 


Decreasing Q 


11 


Benchmark4 


5 


1 


0.5 


2 


Decreasing Q 


y 


Benchmark5 


5 


1 


-1 


2 


Flat Q 


y 


Bcnchmark6 


6 


1 


-1 


2 


Flat Q 


n 


Benchmark? 


5 


1 


-1 


1 


Flat Q 


y 



Table 1. Summary of the benchmarking simulations described in Section[4] p and q are the initial surface mas s density and tem perature 
profiles, E oc R~ p and T oc R~ q , respectively. Simulations Benchmarkl-4 have been set up in the same way as |Rice et al.| ( |2005| whereas 
simulations Bcnchmark5-7 have been set up with a uniform Toomre stability profile over the entire disc. 



density and temperature profiles of the disc are E oc R 1 
and T oc R~ 5 , respectively. The magnitudes of these are set 
such that the Toomre stability parameter (equation [T| at 
the outer edge of the disc, Q m i n = 2. This gives an aspect 
ratio, H/R w 0.05. We model the 1 Mq star in the centre of 
the disc using a sink particle ( Bate et al.|1995 1. At the inner 
disc boundary, particles are accreted onto the star if they 
move within a radius of 0.025 AU of the star or if they move 
into 0.025 ^ R < 0.25AU and are gravitationally bound to 
the star. At the outer edge, the disc is free to expand. 

The simulation was run using a ratio of specific heats, 
7 = 5/3 and hence according to Rice et al. ( 2005[ ), /3 C rit ~ 
6—7. We find that the critical value is « 5.6 since this is 
the lowest value of /3 that the discs can have without frag- 
menting (compare simulations Benchmarkl-3). According to 
equation [4] this is equivalent to a critical value of the gravi- 
tational stress, a?Gi,max ~ 0.07 which is similar to the value 
of ~ 0.06 obtained by Rice et al. (20051. Given the differ- 



ences between the codes, we consider this level of agreement 
acceptable. We therefore compare our remaining simulations 
to this value of /3 or it . 

Figure [l] shows the initial Toomre stability profile of the 



[Tlsho 

(2005 



Rice et al. pOOo] disc (solid line) that is replicated here. As 



a simulation is started, the disc heats up due to the heat- 
ing from gravitational instability, the resulting compression 
and viscous heating, and cools on the cooling timescale de- 
fined by the cooling parameter, j3. Consequently, it is ex- 
pected that the initial disc temperature profile would not 
play a part in the resulting evolution of the disc. We there- 
fore test this by setting up a disc with the same surface 
mass density profile, £ oc R~ , but with a temperature pro- 
file, T oc R, so that its initial Toomre stability profile is 
flat (i.e. constant over the entire disc) with Q = 2 (Fig- 
ure [l] dotted line). These discs were run for /3 = 5 and 6. 
Figure [2] shows the images of the evolved disc with decreas- 
ing Toomre stability profile and the flat Q disc, run with 
a cooling time, p = 5 (simulations Benchmark4 and Bench- 
marks, respectively). The first fragments form at w 20 AU in 
both discs irrespective of the initial temperature profile. Fig- 
ure [3] shows the final Toomre stability profiles of both discs 
run with /? = 6 (simulations Benchmarkl and Benchmark6). 
Neither of these discs fragment and both discs evolve into 
a steady-state with very similar Toomre stability profiles. It 
can be seen that the change in initial temperature profile 
does not make a difference to the final results since with 
/3 — 5 both discs fragment at the same radius and in the 




10 15 
Radius [AU] 

Figure 1. Azimuthally averaged values of the Toomre parameter 
for the initial discs with decreasing Toomre stability profile (sim- 
ulations Benchmarkl-4), set up in the same way as |Rice et al.| 
l |2005| solid line), and with a flat Q profile with Q = 2 (simula- 
tions Benchmark5-7; dotted line). The critical value of Q cr it = 1 
is also marked. 



non-fragmenting cases, the final Toomre stability profiles are 
very similar. Since the temperature in a disc evolves, it is re- 
assuring that the initial temperature profile of the disc does 
not play a part in the outcome. 

As mentioned in Section [I] current wisdom is that ac- 
cording to fragmentation theory, if the Toomre stability pa- 
rameter is below unity and the timescale on which the disc 
cools is faster than a critical value, then the disc should 
fragment. Therefore, if a disc was set up so that its initial 
Toomre stability profile was flat with Q = 1, it would be 
expected that fragments would form everywhere in the disc 
soon after the simulation is started. Figure [4] shows the re- 
sults of this simulation (Benchmark7) . It can be seen that 
despite starting the simulation in a marginally stable state 
where any part of the disc may fragment soon after the simu- 
lation is started, the disc only fragments in the outer regions. 
This illustrates that the disc fragmenting in the outer regions 
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x [AU] x [AU] 

Figure 2. Surface mass density rendered image of the first fragments forming in the simulations using a decreasing Toomre stability 
profile (simulation Benchmark4, left image) and a disc set up with a fiat Q profile with Q = 2 (simulation Benchmark5, right image). 
The discs were run with p = 5. In both cases the discs first fragment at R{ Ri 20 AU confirming that the initial temperature profile 
does not play a part in the evolution of the discs. The colour scale is a logarithmic scale ranging from log E = — 7 (dark) to —3 (light) 
Mq/AU 2 . 




10 15 
Radius [AU] 

Figure 3. Azimuthally averaged values of the Toomre parameter 
for the discs with initially decreasing (solid line) and flat (dotted 
line) Toomre stability profiles (simulations Benchmark!, and 6, 
respectively). The discs were run with /3 = 6. Despite having 
different initial temperature profiles, both discs reach a steady- 
state with very similar Toomre stability profiles, confirming that 
the initial temperature does not play a part in the evolution of 
the discs. The critical value of Q cr it = 1 is also marked. 



< 




Figure 4. Surface density rendered image of the fragmenting 
disc in simulation Benchmark? with an initial flat Q profile with 
Q = 1. The simulation was run with j3 = 5. Despite the initial 
disc being in a state of marginal stability such that, in theory, 
any part of the disc may fragment, the disc only fragments in the 
outer regions. The colour scale is a logarithmic scale ranging from 
log £ = -8 (dark) to -3 (light) M /AU 2 . 



cannot be related to the initial value of the Toomre stability 
profile, Q, and more importantly, fragmentation cannot be 
a function of /3 alone. 
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5 MAIN SIMULATIONS 

In this section, we describe the initial conditions for all the 
individual numerical simulations we have performed to test 
our analytical predictions from Section [2] Table [2] provides 
a summary of the initial conditions as well as the radius at 
which the first fragment forms in the discs that do fragment. 

We set up a series of Reference discs with A/disc 
:'), 1 VI . consisting of 250,000 SPH particles, spanning 0.25 ^ 
R ^ 25AU, surrounding a IMq star, modelled using a sink 
particle. The inner and outer radial disc boundaries have 
been set up in the same way as the benchmarking discs in 
Section [4] All the discs in this section have been set up with 
a flat Q profile. Therefore, as the surface mass density pro- 
file is varied, the initial temperature profile, q, is also varied 
accordingly, where T oc R~ q . The Reference discs have been 
set up so that E oc and T oc R, normalised so that 

Q = 2. We highlight where we have differed from these ini- 
tial conditions in the remaining simulations. 

The analytical work presented in Section[2]suggests that 
for shallow surface mass density profiles, p < 2, fragments 
would form in the outer regions of the discs, whilst for discs 
with steeper surface mass density profiles, p > 2, the disc 
would fragment in the inner regions. We therefore test dif- 
ferent values of the slope of the surface mass density profiles 
using p = 1 (simulations Reference-beta5.5 and Reference- 
beta6), 1.5 (simulations pl.5-beta3.5 and pl.5-beta4), 2.0 
(simulations p2-beta2, p2-beta3 and p2-beta3.5) and 2.5 
(simulations p2.5-beta2, p2.5-beta3.5 and p2.5-beta4). In 
addition, we also carry out a further simulation which is 
the same as simulation p2.5-beta3.5 but with an initial flat 
Q profile, Q = 5 (i.e. so that the initial temperature is 25/4 
times hotter than the disc in simulation p2.5-beta3.5), to 
test the effects of an initially hotter disc on the location of 
fragmentation (simulation p2.5-beta3.5-Q5). 

The analysis also suggests that for a disc with a fast 
enough cooling timescale such that it would fragment, the 
location at which the first fragment would form would move 
inwards to smaller radii as the cooling timescale is decreased. 
We therefore test the effect of decreasing /3 on the fragment 
location by simulating the Reference disc (i.e. with a sur- 
face mass density profile, p = 1) with different values of 
the cooling timescale, /3 = 5.5, 5, 4, 3, 2 and 1 (simula- 
tions Reference-beta5.5, Reference-beta5, Reference-beta4, 
Reference-beta3, Reference-beta2 and Reference-betal, re- 
spectively) . 

We also argued that varying the disc or star mass would 
affect fragmentation. We test the effects of doubling and 
halving the star mass in simulations pl-Mstar2 and pl- 
Mstar0.5, respectively, and compare these to the Reference- 
beta5 simulation. We also carry out extensive tests of the 
effects of varying the disc mass firstly by doubling and 
halving the disc mass (simulations pl-Mdisc0.2 and pl- 
Mdisc0.05, respectively) and secondly by considering more 
extreme disc masses of O.OIMq (simulations pl-beta0.3- 
MdiscO.01, pl-betal-Mdisc0.01, pl-beta2-Mdisc0.01, pl- 
beta2.5-Mdisc0.01 and pl-beta3-Mdisc0.01), O.3M (sim- 
ulation pl-beta8-Mdisc0.3), O.5M (simulation pl-betalO- 
Mdisc0.5) and IMq (simulations pl-beta5-Mdiscl, pl- 
betalO-Mdiscl and pl-betal5-Mdiscl) whilst maintaining a 
central star mass of IMq. 

The analytical work presented in Section [2] also showed 
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Radius [AU] 

Figure 5. Initial surface mass density profiles of the discs used in 
simulations pl-beta6 and pl-beta6-extended. The extended disc 
has the same surface mass density profile as the smaller disc. 



that for a shallow surface mass density profile (p < 2), the 
radius of the disc may be the limiting aspect that causes a 
disc not to fragment. We therefore test a series of extended 
discs which have outer radii, R ou t = 50AU. Simulations pl- 
beta6-extended, pl-beta7-extended and pl-beta8-extended 
are set up so that S and p are the same as in the Reference 
discs (Figure [5|. However, to extend the disc to Rout = 
50AU, the disc masses are increased to Mdisc = 0.2Mq. We 
run this simulation for /3 = 6, 7 and 8 (i.e. values that are 
larger than the critical values identified in Section [IJ. (Note 
that in order to keep the mass of the individual SPH particles 
the same as in the Reference simulations, we use 500,000 
particles to make up this disc). In addition we also set up an 
extended disc with a surface mass density profile, p = 1.5, 
which has a disc mass of 1.5Mq (so that E is the same 
as in simulation pl.5-beta4). We run this simulation using 
/3 = 4 (simulation pl.5-beta4-extended). (As before, since 
we wish to keep the mass of the SPH particles the same as in 
simulation pl.5-beta4, we use 375,000 particles in this disc.) 

Furthermore, we progress the analysis of extended discs 
by simulating two further discs (using 500,000 particles): 
the first is the same as that in pl-beta6-extended but using 
a total disc mass of Mdisc = O.IMq (simulation pl-beta6- 
MdiscO.l-extended) so that the total disc mass is the same 
as in pl-beta6 but E is smaller. The second is also the 
same as in pl-beta6-extended but the central star mass is 
also Mi, — 2Mq so that the disc to star mass ratio is kept 
constant (simulation pl-beta6-Mdisc0.2-Mstar2-extended). 
Both of these discs are run with j3 = 6. 



6 RESULTS 

For the analysis that follows, the key aspect about the frag- 
ments that will be considered will be the first fragment that 
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Simulation 


P 


V 


s 








A^disc 




M diso /M* 


Q 


Disc 


Rf 








name 






[M /(AU) 2 ] 










radius [AU] 


[AU] 








Reference-beta6 


6 


1 


6.4 


X 


lo- 


4 


0.1 


1 


0.1 


2 


25 


- 


- 






Reference-beta5 . 5 


5.5 


1 


6.4 


X 


ur 


-4 


0.1 


1 


0.1 


2 


25 


22 


1.4 x 


1(V 


-2 


Reference-beta5 


5 


1 


6.4 


X 


lo- 


-4 


0.1 


1 


0.1 


2 


25 


20 


1.3 x 


10" 


-2 


Reference-beta4 


4 


1 


6.4 


X 


ur 


-4 


0.1 


1 


0.1 


2 


25 


20 


1.3 x 


ur 


-2 


Reference-beta3 


3 


1 


6.4 


X 


io- 


-4 


0.1 


1 


0.1 


2 


25 


8 


5.1 x 


i(r 


-3 


Reference-beta2 


2 


1 


6.4 


X 


lo- 


4 


0.1 


1 


0.1 


2 


25 


3 


1.9 x 


io- 


-3 


Reference-betal 


1 


1 


6.4 


x 


ur 


-4 


0.1 


1 


0.1 


2 


25 


2.5 


1.6 x 


10- 


-3 


pl.5-beta4 


4 


1.5 


1.8 


X 


10" 


-3 


0.1 


1 


0.1 


2 


25 


_ 


_ 






pl.5-beta3.5 


3.5 


1.5 


1.8 


x 


ur 


-3 


0.1 


1 


0.1 


2 


25 


18 


7.5 x 


iir 


-3 


pl.5-beta3 


3 


1.5 


1.8 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


1.7 


2.3 x 


iir 


-3 


p2-bota3.5 


-i 5 


2 


3.5 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 










p2-beta3 


3 


2 


3.5 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


0.45 


3.5 x 


io- 


-3 


p2-beta2 


2 


2 


3.5 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


0.3 


3.5 x 


iir 


-3 


p2.5-beta4 


4 


2.5 


1.4 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


_ 


_ 






p2.5-beta3.5 


3.5 


2.5 


1.4 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


0.4 


7.0 x 


io- 


-3 


p2.5-beta3.5-Q5 


3.5 


2.5 


1.4 


x 


10" 


-3 


0.1 


1 


0.1 


5 


25 


0.3 


8.1 x 


io- 


-3 


p2.5-beta3 


3 


2.5 


4.4 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


0.3 


8.1 x 


iir 


-3 


p2.5-beta2 


2 


2.5 


4.4 


x 


10" 


-3 


0.1 


1 


0.1 


2 


25 


0.35 


7.5 x 


10- 


-3 


pl-Mstar2 


5 


1 


6.4 


X 


10" 


4 


0.1 


2 


0.05 


2 


25 


- 


- 






pl-Mstar0.5 


5 


1 


6.1 


X 


10" 


-4 


0.1 


0.5 


0.2 


2 


25 


13 


1.7 x 


io- 


-2 




5 


\ 


1.3 


X 


10" 


3 


2 




2 


2 




14 


1.8 x 


10" 


-2 


pl-Mdisc0.05 


5 


1 


3.2 


X 


ur 


-4 


0.05 


1 


0.05 


2 


25 










pl-betal-Mdisc0.01 


1 


1 


6.4 


X 


io- 


-5 


0.01 


1 


0.01 


2 


25 


6.5 


4.2 x 


io- 


-4 


pl-beta2-Mdisc0.01 


2 


1 


6.4 


X 


ur 


-5 


0.01 


1 


0.01 


2 


25 


15 


9.6 x 


io- 


-4 


pl-beta2.5-Mdisc0.01 


2.5 


1 


6.4 


X 


ur 


-5 


0.01 


1 


0.01 


2 


25 


17 


1.1 X 


iir 


-3 


pl-beta3-Mdisc0.01 


3 


1 


6.4 


X 


10" 




0.01 


1 


0.01 


2 


25 


- 


- 






pl-beta8-Mdisc0.3 


8 


1 


1.9 


X 


io- 


-3 


0.3 


1 


0.3 


2 


25 


- 


- 






pl-betalO-Mdisc0.5 


10 


1 


3.2 


X 


io- 


3 


0.5 


1 


0.5 


2 


25 










pl-beta5-Mdiscl 


5 


1 


6.4 


X 


iir 


-3 


1 


1 


1 


2 


25 


5.5 


3.5 x 


10" 


-2 


pl-beta7-Mdiscl 


7 


1 


6.4 


X 


ur 


-3 


1 


1 


1 


2 


25 










pl-botalO-Mdiscl 


10 


1 


6.4 


X 


io- 


-3 


1 


1 


1 


2 


25 










pl-betal5-Mdiscl 


15 


1 


6.4 


X 


io- 


-3 


1 


1 


1 


2 


25 










pl-beta6-extended 


6 


1 


6.4 


X 


io- 


4 


0.2 


1 


0.2 


2 


50 


24.5 


1.6 x 


io- 


-2 


pl-beta7-extended 


7 


1 


6.4 


X 


io- 


-4 


0.2 


1 


0.2 


2 


50 


29 


1.9 x 


io- 


-2 


pl-beta8-extended 


8 


1 


6.4 


X 


iir 


-4 


0.2 


1 


0.2 


2 


50 


30 


1.9 x 


iir 


-2 


pl.5-beta4-extended 


4 


1.5 


1.8 


X 


ur 


-3 


0.15 


1 


0.15 


2 


50 


33 


1.0 x 


iir 


-2 


pl-beta6-MdiscO. 1-extended 


6 


1 


3.2 


X 


io- 


-4 


0.1 


1 


0.1 


2 


50 


10 


1.3 x 


10" 


-2 


pl-beta6-Mdisc0.2-Mstar2-extended 


6 


1 


6.4 


X 


io- 


4 


0.2 


2 


0.1 


2 


50 


31 


1.1 x 


io- 


-2 



Table 2. Summary of the main simulations, p describes the initial surface mass density profile, £ oc R~ p , and S is the normalisation 
constant required to produce a disc with mass M^ ac . The final column represents the RHS of equation [7] for the location at which the 
first fragment forms, R{. The simulations were run with an initial flat Toomre stability profile, Q. 



forms. This is because subsequent evolution of the disc fol- 
lowing an initial fragmentation stage may involve additional 
complexities that are beyond the scope of this paper. Ta- 
ble [2] summarises the key fragmenting results. The radius at 
which the first fragments form, Rf, has been determined by 
eye from the disc images. It is important to note that as seen 
in past simulations (e.g. Lodato & Rice 2004, Meru & Bate 



20101, the surface mass density profile does not change sig- 
nificantly during the simulations, particularly for low mass 
discs (Mdisc < 0.2M©). We highlight where the surface mass 
density profiles do change significantly and discuss the ef- 
fects of this. 



6.1 Fragmentation dependency on the surface 
mass density profile 



Figures |6|9| show that as the surface mass density profile 
steepens, the location at which the first fragment forms 
moves to smaller radii in the disc. The analytical theory 
presented in Section [2] shows that for a shallow surface mass 
density fall off where p < 2, the fragments form in the outer 
regions of the disc (provided the cooling criterion is also sat- 
isfied) . This is indeed the case for simulations with p = 1 and 
p = 1.5 (simulations Reference-beta5.5 and pl.5-beta3.5, re- 
spectively) as Figures |6] and [T] show that the fragments form 
at Rf w 20 AU and w 19 AU, respectively. Figures [To] and [TT| 
show the radial profile of the aspect ratios (calculated by az- 
imuthally averaging the sound speed at each radii and divid- 
ing by QR at that radii) in the discs compared to the RHS 
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-20 20 

x [AU] 

Figure 6. Surface mass density rendered image of the fragment- 
ing disc with initial surface mass density profile S oc R~ The 
simulation (Reference-beta5.5) used (3 = 5.5. The fragment forms 
in the outer regions of the disc, confirming the analytical predic- 
tions in Section [2] The colour scale is a logarithmic scale ranging 
from log S = -6 (dark) to -3 (light) Mq/AU 2 . 




-20 20 

x [AU] 

Figure 7. Surface mass density rendered image of the fragment- 
ing disc with initial surface mass density profile S oc i? -3 / 2 . The 
simulation used j3 = 3.5. The fragment forms in the outer regions 
of the disc, confirming the analytical predictions in Section[2] The 
colour scale is a logarithmic scale ranging from log E = — 7 (dark) 
to -2 (light) M Q /AU 2 . 

of equation [7] at a time shortly before the discs fragment. 
It can be seen that condition [7] is satisfied at the region in 
which the first fragment forms shortly after. The oscillations 
in H/R are due to temperature fluctuations since although 
the cooling rate in the disc changes smoothly with radius, 
the heating of the disc occurs primarily in the spiral shocks. 
This therefore confirms the analytical predictions for shallow 



surface mass density profiles (p < 2) presented in Section [2] 
It is important to note that for a flat Q profile, the tem- 
perature profiles in the discs are an increasing function of 
radius (T oc R) for p = 1 and a constant temperature profile 
for p = 1.5, yet the discs still fragment in the outer regions 
(again, re-emphasising that the initial temperature profile 
does not play a part in the disc evolution). 

The analytical theory for p > 2 suggested that if the disc 
was to fragment, it would do so in the inner regions of the 
disc. Figures [8] and [9] show that this is indeed the case. Fig- 
ure |12| shows that the analytical condition is just satisfied for 
simulation p2-beta3 at the location at which the fragment 
forms. Figure fl3] shows that for simulation p2.5-beta3.5, the 
analytical condition is also satisfied at the location at which 
fragmentation occurs soon after. 

We therefore show that as the surface mass density pro- 
file is steepened so that more of the mass is concentrated 
in the inner regions of the disc, fragmentation moves to- 
wards smaller radii. It is important to note that the trend 
that fragmentation moves to smaller radii for steeper surface 
mass density profiles is valid even when considering a uni- 
form value of P (compare simulations Reference-beta3, pi. 5- 
beta3, p2-beta3 and p2.5-beta3 which are run using /3 = 3 
and fragment at ~ 8, 1.7, 0.45 and 0.3 AU, respectively). 

In addition, the results summarised in Table [2] show 
that a single value of /3 cr i t is not applicable over all surface 
mass density profiles since the minimum value of ft that a 
disc can have without fragmenting varies with the surface 
mass density profile. 

Simulation p2.5-beta3.5-Q5 was the same as simulation 
p2.5-beta3.5 but had an initial disc that was hotter by a 
factor of 25/4. The results of this simulation show that the 
disc still fragmented in the inner regions. 

6.2 Effect of the cooling timescale, /3, on the 
fragment location 

In Section [2] the analytical work presented suggested that 
for a fragmenting disc with a shallow surface mass density 
profile (p < 2), a decrease in the value of f3 would cause the 
location at which the first fragment forms to move inwards 
to smaller radii. 

Figure [14] shows the radius at which the first frag- 
ment forms for different values of j3 (simulations Reference- 
beta5.5, Reference-beta5, Reference-beta4, Reference-beta3, 
Reference-beta2 and Reference-betal). We can see a clear 
trend showing that the radius of fragmentation moves in- 
wards for more efficient cooling. 

6.3 The influence of star mass on fragmentation 

In Section |2j we showed that decreasing the star mass is 
more likely to cause conditons[6]and[7]to be satisfied over a 
larger part of the disc and hence the disc is more likely to 
fragment. We test three identical discs with star masses of 
0.5, 1 and 2 Mq (simulations pl-Mstar0.5, Reference-beta5 
and pl-Mstar2, respectively) which are run with the same 
cooling timescale, /3 = 5. It can be seen from Table[2]that the 
Reference-beta5 disc first fragments at R{ « 20 AU. How- 
ever, when the star mass is halved, the first fragment forms 
at Rf » 13 AU. Since the RHS of condition [7] oc R 2 /M ir , if 
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Figure 8. Surface mass density rendered image of the fragmenting disc with initial surface mass density profile E oc R~ 2 (simulation 
p2-beta3). The simulation used /3 = 3. The fragment forms in the inner regions of the disc as shown by the zoomed in image of the disc, 
confirming the analytical predictions in Section [2] The colour scale is a logarithmic scale ranging from log E = — 11 (dark) to 2 (light) 
Mq/AU 2 in the zoomed out image and from log E = —3.5 (dark) to 1 (light) Mq/AU 2 in the zoomed in image. 



20 



i 



-20 




Figure 9. Surface mass density rendered image of the fragmenting disc with initial surface mass density profile E oc i? -5 / 2 (simulation 
p2.5-beta3.5). The simulation used ft = 3.5. The fragment forms in the inner regions of the disc, confirming the analytical predictions in 
Section|2] The colour scale is a logarithmic scale ranging from log E = —12 (dark) to 3 (light) Mq/AU 2 in the zoomed out image and 
from log E = —4 (dark) to 0.4 (light) Mq/AU 2 in the zoomed in image. 



the star mass is halved for the same value of j3 (and hence 
the same value of the LHS of condition [Tf, the radius at 
which the first fragment forms, Rf, decreases by a factor of 
y/2. Conversely, doubling the star mass makes it harder for 
the condition to be satisfied and indeed the disc in simula- 
tion pl-Mstar2 does not fragment. Instead it settles into a 
state of marginal stability with Q « 1. 



6.4 The influence of disc mass on fragmentation 

The analytics presented in Section [2] showed that increas- 
ing the disc mass (and hence increasing E ) allows condi- 
tons[6]and[7]to be satisfied over a larger part of the disc and 
hence the disc is more likely to fragment. We initially test 
this by comparing the results of simulations pl-Mdisc0.05, 
Reference-beta5 and pl-Mdisc0.2 which are identical discs 
except that the disc masses are 0.05Mq, O.IMq and 0.2Mq, 
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Figure 10. Plot of disc aspect ratio, H/R, (solid line), against the 
RHS of equation]?] (dotted line) for simulation Reference-beta5.5. 
Condition[7]is satisfied at « 20AU where the disc first fragments, 
confirming the analytical predictions in Section [2] 
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Figure 11. Plot of disc aspect ratio, H/R (solid line), against 
the RHS of equation [7] (dotted line) for simulation pl.5-beta3.5. 
Condition[7]is satisfied at « 19AU where the disc first fragments, 
confirming the analytical predictions in Section [2] 



respectively. Table[2]shows that doubling the disc mass from 
O.IM0 to 0.2Mq does indeed cause the fragmentation condi- 
tions derived in Section[2]to be satisfied over a larger portion 
of the disc since the first fragments form at Rf w 20 AU and 
~ 14 AU, respectively. Since the RHS of condition [7] oc T,R 2 , 
if the disc mass (and hence E) is doubled for the same value 



0.02 



0.01 




10 15 
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0.01 - 




0.6 
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Figure 12. Plot of disc aspect ratio, H/R (solid line), against the 
RHS of equation [7] (dotted line) for simulation p2-beta3 for the 
radial range of the entire disc (upper panel) as well as zoomed 
into the inner regions (lower panel). Condition [7] is marginally 
satisfied at Ri 0.4AU where the disc first fragments, confirming 
the analytical predictions in Section [2] 



of f3 (and hence the same value of the LHS of condition [7]) , 
the radius at which the first fragment forms, Rf, decreases 
by a factor of v2. However, halving the disc mass makes it 
harder for the conditions to be satisfied and consequently, 
the disc does not fragment. 

In addition, we also simulate a very low mass disc 
(Afdisc = O.OIM0) and found that it fragments if f3 = 
1, 2 and 2.5 but not for /3 = 3 (simulations pl-betal- 
MdiscO.01, pl-beta2-Mdisc0.01, pl-beta2.5-Mdisc0.01 and 
pl-beta3-Mdisc0.01, respectively). As /3 increases, the frag- 
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Figure 13. Plot of disc aspect ratio, H/R, (solid line), against the 
RHS of equation[7] (dotted line) for simulation p2.5-beta3.5 for the 
radial range of the entire disc (upper panel) as well as zoomed 
into the inner regions (lower panel) . Condition [7] is satisfied at 
Rj 0.4AU where the disc first fragments, confirming the analytical 
predictions in Section [2] 




Figure 14. The radius at which the first fragment forms in the 
Reference simulations. The discs in these simulations are identical 
with a surface mass density profile, p = 1, but were run with 
different values of the cooling timescale in units of the orbital 
timescale, /3. The radius at which the first fragment forms moves 
inwards with more efficient cooling. 




0.5 

log Radius [AU] 



ment location moves out in the disc, as found in Section [6. 2| 
It is clear that, as with varying the star mass, the disc mass 
plays a crucial role in the fragmentation and the condition 
for fragmentation cannot simply be described using a single 
critical value of the cooling timescale. 

In addition, we simulate higher mass discs, Mdi BC = 
O.3M0 and 0.5, which are run using /3 = 8 and 10 re- 
spectively (simulations pl-beta8-Mdisc0.3 and pl-betalO- 
Mdisc0.5) as well as discs with Mdi sc = IMq which we sim- 
ulate using f3 = 5, 10 and 15 (simulations pl-beta5-Mdiscl, 



Figure 15. Surface mass density profiles for simulation pl-beta7- 
Mdiscl at the start (solid line) and at a time more than 4 ORPs 
later (dotted line). Unlike the low mass simulations whose surface 
mass density profiles do not change throughout the simulations, 
the profile for this disc steepens causing a change in the effective 
values of S and p. 
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Figure 16. Plot of disc aspect ratio, H/R (solid line), for sim- 
ulation pl-beta7-Mdiscl, against the RHS of equation [7] using 
the initial values of S D and p (dashed line) and the new values 
of S c and p determined after the disc has evolved for > 40RPs 
by which time its surface mass density profile has changed. The 
condition is satisfied using the initial values of E D and p but not 
using the new values and hence the disc does not fragment. 




-40 -20 20 40 



x [AU] 

Figure 17. Surface mass density rendered image of the frag- 
menting disc in simulation pl-beta8-extended with initial surface 
mass density profile £ <x R~ , but extending to 50 AU rather 
than 25 AU. This simulation was run with = 8. According to 
|Rice et al.| ( |2005[ ) , this disc should not fragment since the cooling 
timescale f} is larger than the critical value previously obtained 
with a radius of 25 AU. This simulation shows that the frag- 
mentation criterion is more complex than a single critical cooling 
parameter. The colour scale is a logarithmic scale ranging from 
log S = -8 (dark) to -2 (light) M /AU 2 . 



pl-betalO-Mdiscl and pl-betal5-Mdiscl, respectively). We 
find that with the exception of simulation pl-beta5-Mdiscl, 
the discs do not fragment. Figure [15] shows the surface mass 
density profile of simulation pl-beta7-Mdiscl at the start 
and more than 4 ORPs after the start of the simulation. 
It can be seen that unlike the lower mass discs, the profile 
steepens and the value of E increases. This is the case for all 
the non-fragmenting high mass disc simulations. Figure [16] 
shows the plot of the aspect ratio profile of this simulation 
against the RHS of condition [7] and shows that condition [7] is 
just satisfied in the inner regions. However, since during the 
simulation the disc mass redistributes itself, the surface mass 
density profile changes and consequently, using the newly 
obtained values of E and p, condition [7] is not satisfied. We 
note that the only high mass disc that does fragment (simu- 
lation pl-beta5-Mdiscl), does so because the cooling time is 
so rapid that fragmentation occurs before the disc has had 
the chance to restructure itself. 



6.5 The role of the disc radius on fragmentation 

In Section[2] we showed that for shallow surface mass density 
profiles (p < 2), fragmentation might occur for any value of 
the cooling timescale, if the disc is large enough. 

Simulation Reference-beta6 (a disc with R out = 25 AU) 
did not fragment and though we did not run the same sim- 
ulation with ft — 7 or 8, we would expect that they would 
also not fragment. However, extended discs with the same 
values of E and p as Reference-beta6 do indeed fragment 
for j3 — 6 (simulation pl-beta6-extended), f3 = 7 (pl-beta7- 
extended) and ft = 8 (pl-beta8-extended; Figure 17 1 with 



the first fragments forming at Rf « 25, 29 and 30 AU, re- 
spectively. Similarly, we also simulate an extended disc with 
p — 1.5 and ft — 4 (simulation pl.5-beta4-extended) and 
show that whilst the same disc truncated at R out = 25AU 
does not fragment, the extended disc does indeed fragment 
(at R { « 33 AU). 

In addition, given that in Section |6.4| we showed that 
the disc mass plays a part in whether fragmentation occurs 
or not, we simulate a O.IMq disc which extends to R ou t = 
50AU (simulation pl-beta6-Mdisc0.1-extended). The sur- 
face mass density profile is the same as in simulation pl- 
beta6-extended (p — 1) but E is decreased. The results 
show that the disc fragments at R ~ 40AU (c.f. ~ 25 AU 
for simulation pl-beta6-extended). Therefore, whilst the disc 
mass affects where in the disc the first fragment forms, 
the conclusion that a small disc which does not fragment 
for a particular value of /3 may fragment at larger radii 
for the same value of /3 is still valid. Furthermore, we can 
see that if the disc to star mass ratio is kept constant at 
M^isc/Mi, = 0.1 for an extended disc, the disc fragments 
(at Rf w 34 AU; simulation pl-beta6-Mdisc0.2-Mstar2- 
extended) whilst the small disc with the same disc to star 
mass ratio does not fragment, further corroborates the fact 
that the radius of the disc is important. 



7 DISCUSSION 

It has previously been accepted that for a self-gravitating 
disc whose only source of heating is internally generated 
from the gravitational instability, the disc will fragment if 
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the cooling timescale is short enough ( Gammie|2001 1 . How- 
ever, we find that fragmentation at a given radius is not 
only dependent on the cooling timescale, ft, but also on the 
disc surface density (i.e. disc mass and profile) and the star 



This is in contrast to Rice et al. ( 2005 1 who suggested 



that the fragmentation criterion is independent of the disc 



mass though in agreement with Rice et al. ( 2003 1 who found 



that for a higher disc mass, fragmentation was easier: using 
/3 = 5, they found that a O.25M0 disc fragmented whilst a 
O.IM0 disc did not, but instead required a lower value of 
ft. In particular, we highlight that in the past, it has been 
thought that a massive disc is required for fragmentation to 
occur. However, we show that it is indeed possible for low 
mass discs (Mdisc ~ O(O.OI)Mq) to fragment if the cooling 
in the discs is rapid enough. On the other hand, for high 
mass discs (A/disc ^ O.3M0 within R ou t = 25 AU), the discs 
do not fragment, unless the cooling time is fast, due to a 
steepening of the surface mass density profile and an in- 
crease in E D making condition [7] harder to satisfy. |Lodato| 
& Rice ( 2005 ) also find that with larger disc masses, matter 
is redistributed causing the surface mass density profile to 
steepen. 

In addition, we find that the critical value of ft found 
for one particular surface mass density profile is not ap- 
plicable to another disc mass distribution. Our simulations 
show that for a steeper surface mass density profile, the 
cooling timescale req uired for a disc to fragment is smaller. 



Cossins et al. 



Pcrit 
/?crit 



6 



|2010|) found that for a dis c with E oc R~ 3/2 , 
( 2005) found that 



Rice et al. 



— 4.5 whereas 
7 for a disc with E oc 



7? . In addition, 



Rice 



et al. ( 2003 1 found that for a surface mass density profile of 
E oc R~ ' /4 , the fragmentation boundary for a 0.1M Q disc 
around a 1M© star is between ft = 3 and 5 which [Rice et al."| 
( 2005 1 note is inconsistent with their results. However, our 



results explain these previous inconsistencies present in the 
literature. 

We also find that for p < 2, if a disc does not frag- 
ment for a particular cooling timescale in units of the 
orbital timescale, ft, a larger disc with the same surface 
mass density profile may well fragment (compare simula- 
tions Reference-beta6 and pl-beta6-extended or pl.5-beta4 
and pl.5-beta4-extended, and also pl-beta7-extended or pl- 
beta8-extended) . Therefore, a critical cooling timescale can 
only be specified for a particular surface mass density at 
a particular disc radius and can therefore not be a general 
rule. The previous fragmentation criterion found by |Rice| 
|et al.| ( |2005[ ) was for a disc to star mass ratio of 0.1 with 
M* = IMq, a surface mass density profile of E oc -R -1 and 
R ou t = 25 AU. 

Finally, we also show that the mass of the star plays a 
part in whether a disc fragments or not. Therefore, as shown 
by the RHS of condition [7j whether a disc fragments or not 
is essentially a "trade-off" between the local surface mass 
density and the star mass. 

7.1 The link between ft, Mdi sc , M+, the surface 
mass density profile, p, and fragmentation 

Section [6] shows that the fragmentation criterion is clearly 
a complex problem which cannot simply depend on a single 
critical cooling timescale as has previously been thought to 




Figure 18. Logarithmic graph showing the trend between ft and 
ERj /M* determined by considering the location at which the first 
fragment forms in the discs, Rf. The results include those simula- 
tions with a surface mass density profile, p = 1 (filled triangles), 
p = 1.5 (open triangles), p = 2 (open squares) and p = 2.5 
(crosses). It is clear that a single critical value of ft is not the case 
for all discs and that there is a relation between ft, Mdi sc , M* 
and the surface mass density profile, p, that determines whether 
fragmentation occurs or not. The trendline has been determined 
by considering discs with shallow surface mass density profiles, 
p < 2 only as those discs with p > 2 will always fragment in the 
innermost regions first. The grey shaded region is where we expect 
subsequent fragmentation may take place in discs with p < 2. 



be the case. Equation[7]and the results presented here clearly 
show that there is a link between the cooling timescale (in 
terms of /3), the disc mass (or more accurately, the local 
surface mass density) and the star mass. 

Such a link can be explained physically. As a disc cools, 
gravitational instability develops resulting in density fluc- 
tuations above and below the unperturbed density, 5p/p. 
The spiral structures involve shocks which produce heat that 
may balance the disc's cooling, thus reaching an equilibrium 
state. If the disc mass was irrelevant, for any particular star 
mass, and ft, one would expect the fluctuations, Sp/p, to be 
the same in all discs with the same surface density profile. 
However, comparing a low mass disc with a high mass disc 
with the same relative density fluctuations {Sp/p), the den- 
sity enhancement, Sp, in the higher-mass disc will clearly be 
greater. At some disc mass, this enhancement will be self- 
gravitating (i.e. it will be a fragment), while in a lower-mass 
disc the density fluctuation will not form a fragment (unless 
the value of ft is lowered). Similarly, if the disc is kept the 
same and the star's mass is increased, a given density en- 
hancement may be sheared apart by the differential rotation 
so that a lo w v alue of ft will be required for fragmentation, 
shows a graph of ft against Eiif/M* 
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where 

Rf is the radius at which the first fragment forms, and in- 
cludes all the fragmenting simulations presented in this pa- 
per (including the Benchmarking simulations). When inter- 
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preting these results, one should also note that an increase 
in Ei? 2 /M* does not necessarily imply an increased disc to 
star mass ratio: it is possible to have a low disc to star mass 
ratio but using a steeper surface mass density profile. We 
can see that there is a clear trend that as the RHS of equa- 
tion [7] increases, so too does the value of j3 that will allow a 
fragment to form. The trendline presents the relation: 



shown by Rice et al. ( 2005 1 , the ratio of specific heats plays 



AL 



(9) 



where 8 « 1/2 and the constant of proportionality, r\ w 47, 
which we find using a least squares fit. We have performed a 
few calculations with higher resolution (2 x 10 6 SPH parti- 
cles) and find that this trend of j3 increasing with Ei? 2 /M* is 
maintained, although the exact values of r\ and 5 may change 
slightly if we were able to perform all the calculations with 
high resolution. 

It is also important to note that the trendline has only 
been produced using the results of simulations with p < 2. 
This is because for p > 2, the disc will always fragment in 
the innermost regions first. Consequently, if the results for 
p > 2 were included, this would cause the trendline to be 
somewhat skewed. 

The trendline can give some very useful information for 
those simulations with p < 2. The grey region is where we 
predict subsequent fragmentation is feasible. Traversing the 
plot in a vertical direction downwards from the trendline 
into the grey region, at any particular value of Ei? 2 /M* the 
disc will fragment at all values of ft less than the limit given 
by the trendline (though the radius being considered will 
not necessarily be the first location at which fragmentation 
occurs). Similarly, traversing the plot in a horizontal direc- 
tion from the trendline to the right side of the plot into the 
grey region, for a particular value of /3, fragmentation is pos- 
sible at a particular radius if the disc mass is increased or 
star mass is decreased. Similarly, for a particular value of 
P and any one combination of the disc to star mass ratio, 
fragmentation is also possible at larger radii than the loca- 
tion specified by the trendline i.e. to the right hand side of 
the trendline. Therefore, the trendline predicts the minimum 
possible radius at which fragments could in theory form in 
discs with shallow surface mass density profiles, p < 2. 

For discs with p > 2 which fragment in the inner regions, 
we expect that subsequent fragmentation may take place 
further out in the disc as far out as given by the trendline in 
Figure [18] In other words, for these discs, we expect there to 
be a maximum radius outside of which fragmentation will 
not occur (since the surface mass density fall-off is steep, 
the outer regions of these discs may struggle to have enough 
mass for gravitational instability to be significant). However, 
since we stop the simulations soon after the first fragment 
forms due to the increased computational resources required 
to follow the simulations further, we are unable to test this 
prediction. Further work needs to be done in this area and 
is beyond the scope of this paper. We also note that in real 
discs, for fragmentation to occur in the inner regions, the 
cooling time would have to be very small since the dynamical 
time at small radii would also be very small. Such short 
cooling times may not be possible in real discs. 

It is also important to note that these simulations have 
been carried out using a ratio of specific heats, 7 = 5/3. As 



a key role in the fragmentation boundary. We therefore also 
anticipate a further dependency on the equation of state. 
Furthermore, for high mass discs (Mdi sc > 0.3Mq), we find 
that the initial surface mass density conditions (i.e. E and 
p) cannot be used to determine whether the disc will frag- 
ment or not. This is because as the disc evolves, the surface 
mass density profile steepens causing E and p to change. 
Consequently, though parts of the disc may start off in the 
grey shaded region of Figure [18] and hence may be expected 
to fragment, the disc may restructure itself on a timescale 
faster than the cooling timescale such that it moves out of 
this region and hence does not fragment. 



Following the work of Gammie ( 2001 1 and Rice et al 



1:2005 ) , a number of authors have used the critical cooling 
timescale, /3 cr it ( or gravitational stress, CKGl.max) to predict 



fragmentation in realistic discs (e.g. Clarke 2009 Rafikov 



2009||Cossins et al.|2010||Kratter et al.|2010[ ). In light of the 



new results presented here, we would encourage previous 
conclusions based on these critical values to be revisited. 



8 CONCLUSIONS 

We present an analytical approach to examine the fragmen- 
tation of self-gravitating protoplanetary discs, and confirm 
the results using global three-dimensional numerical simula- 
tions. Our key result is that fragmentation does not simply 
depend on the disc cooling timescale, /3, but also on the ratio 
of the surface mass density at radius, R, to the stellar mass, 
i.e. TjR 2 /M^. We find that fragmentation occurs when 



P < vl 



(10) 



where 5 « 1/2 and 77 « 47. For a power-law surface mass 
density, E oc R~ p , this relation predicts the innermost radii 
at which subsequent fragments can form in a disc with shal- 
low surface mass density profiles (p < 2) as well as the radius 
outside of which fragments cannot form in discs with steep 
surface mass density profiles, (p > 2). Generally, we find 
that an increase in the steepness of the disc surface mass 
density profile promotes fragmentation at smaller radii. 
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